%==============================================================================
% This code computes equilibrium with and without CBDC under different CBDC 
% as the economy gets cashless, i.e., the department stores that accept
% both cash and electronic payment moves online and accept only electronic
% payment
%  COPY RIGHT: Yu Zhu
%==============================================================================



clear;close all
options = optimset('maxiter',5e5,'tolX',1e-10,'tolfun',1e-10,'maxfunevals',5e5);

%==========================================
%Define parameters and functions
%==========================================
name = 'LN87_08_final'
load(['pars_post_' name '.mat'])
location='figures/'

alpha1=pars(4);alpha2=pars(5);alpha3=pars(6);beta=pars(7);epsilon=0.001;sigma=pars(1);
cost=pars(8);                   
B=pars(2);
A=pars(10);  %productivity
eta=pars(9);                                       %bank's cost to manage deposits
theta=pars(11);
Nb=pars(3); %number of banks
rq_ratio=pars(12);
i_reserve = pars(13);
year_option = pars(14); % 1: using banking data from 2007 and 2008. 0: using banking data from 2014-2019

%% Define the name of the saved files: last two number indicates reserve requirement ratio and the interest on reserves
if year_option==1
    pi_m = 0.03346;pi_h=pi_m;
    name = [name,'_', num2str(floor(rq_ratio*100)),'_',num2str(floor(i_reserve*100)) '_07_08'];
else
    pi_m = 0.01515;pi_h=pi_m;
    name = [name,'_', num2str(floor(rq_ratio*100)),'_',num2str(floor(i_reserve*100))];
end



picbdc=pi_m;

pi_eff = (1+pi_m)/(1+i_reserve)-1;

nrhogrid =10000;
i_m = 1/beta*(1+pi_m)-1;i_h = 1/beta*(1+pi_h)-1;
u = @(x) ((x+epsilon).^(1-sigma)-epsilon.^(1-sigma))./(1-sigma); %DM utility
du = @(x)(x+epsilon).^(-sigma); %derivative of DM utility
d2u=@(x)(x+epsilon).^(-sigma-1)*(-sigma);
c=@(x)x;                                                         % DM cost
dc=@(x)1;%   derivative of DM cost
d2c=@(x)0;
y_star = fminsearch(@(x)(du(x)-dc(x)).^2,1);                     % y_star
g_func=@(x) (1-theta).*u(x)+theta.*c(x);                         %trading proctocol 
dg_func=@(x) (1-theta).*du(x)+theta.*dc(x);
d2g_func=@(x)(1-theta).*d2u(x)+theta.*d2c(x);
y_grid = 0:0.001:y_star;
p_grid = g_func(y_grid);
p_star = p_grid(end);
y_func = @(p)interp1(p_grid,y_grid,min(p,max(p_grid)));          %trading quantity 
p_func=@(x)min(x,max(p_grid));                                   %trading payment
lambda_func = @(z) max(du(y_func(z))./dg_func(y_func(z))-1,0);          % liquidity premium

f=@(k)A*k.^(eta);%production
df=@(k)A*eta*k.^(eta-1); %marginal production
df_inv = @(rho)(A*eta./(1+rho)).^(1/(1-eta));


dlambda_func=@(z)((d2u(y_func(z)).*dg_func(y_func(z))-du(y_func(z)).*d2g_func(y_func(z)))./dg_func(y_func(z)).^3).*(z<=p_grid(end));
dinv=@(z)beta*(alpha2*lambda_func(z)+1);
ddinv=@(z)beta*alpha2*dlambda_func(z);
z_grid = linspace(0,p_grid(end),100000);
lambda_grid = lambda_func(z_grid);
gamma_grid=linspace(0,20,nrhogrid);
x0=[0,0,0];
x0_R=x0;
opt=1;
icbdc_grid=linspace(-0.01,0.05,1000);
D_grid = linspace(0,max(p_grid),5e3);
pass_variables
ygrid = D_inv_demand(i_m,D_grid,func_var);
alpha2_grid = alpha2 + linspace(0,0.03,500);
alpha3_grid = alpha3 - linspace(0,0.03,500);
eq_picbdc=(1+picbdc)/(1+0)-1;

for i=1:length(alpha2_grid)
    i
    alpha2=alpha2_grid(i);
    alpha3 = alpha3_grid(i);
    % Compute equilibrium without CBDC
    pass_variables
    if i==1
        y=SS_eq_noCBDC(func_var,pi_eff,Nb,ygrid);
    else
        y=SS_eq_noCBDC_fast(func_var,pi_eff,x0,opt);
    end
    
    i_m=(1+pi_m)/beta-1;
    psi= beta*(alpha2*lambda_func(y(2))+alpha3*lambda_func(y(1)+y(2)))+beta;
    drate(i)=1/psi-1;
    loan= fzero(@(x)df(x)-1-y(end),[1e-10,100]);
    
    eq_final_nocbdc(i,:)=[y,drate(i),loan];
    x0=y;
   
    %===============================================
    %Compute equilibrium with CBDC
    %===============================================
    icbdc_real  = (1+picbdc)/(0+1)/beta-1;
    if icbdc_real>psi/beta-1
        eq_cbdc(i,:)=[y,drate(i),loan];
        drate_cbdc(i)=drate(i);
        eliquidity(i,1) = y(2);
    else
        tic
        z_eq= SS_eq( alpha1,alpha2,alpha3,z_grid,lambda_grid,i_m,icbdc_real);
        
        lsupply_max=(1-rq_ratio)*(1+icbdc_real)*beta*z_eq(1,2);
        rho_hat=max(1/(1+pi_eff)-1, ((0+1)/(1+picbdc)-rq_ratio/(1+pi_eff)+cost)/(1-rq_ratio)-1);
        ld_low= fzero(@(x)df(x)-1-rho_hat,[1e-10,100]);
        drate_cbdc(i)=(0+1)/(1+picbdc)-1;
        toc
        if ld_low<lsupply_max
            deposit =  ld_low/(icbdc_real+1)/beta/(1-rq_ratio);
            eq_cbdc(i,:)=[z_eq(1),deposit,rho_hat,drate_cbdc(i),ld_low];
        else
            rhocbdc = df(lsupply_max)-1;
             
            eq_cbdc(i,:)=[z_eq,rhocbdc,drate_cbdc(i),lsupply_max];
            
        end
        eliquidity(i,1) = z_eq(2);

    end
    
    
    
    
end
Output_noCBDC = alpha1*eq_final_nocbdc(:,1)+alpha2_grid(:).*eq_final_nocbdc(:,2)+alpha3_grid(:).*min(eq_final_nocbdc(:,1)+eq_final_nocbdc(:,2),p_star)...
    +B+eq_final_nocbdc(:,end)+f(eq_final_nocbdc(:,end))-eq_final_nocbdc(:,2)...
    +(eq_final_nocbdc(:,2)./(1+eq_final_nocbdc(:,4))-eq_final_nocbdc(:,end))*(1+i_reserve)/(1+pi_m);
Output_CBDC = alpha1*eq_cbdc(:,1)+alpha2_grid(:).*eliquidity+alpha3_grid(:).*min(eq_cbdc(:,1)+eliquidity,p_star)+B+eq_cbdc(:,end)+f(eq_cbdc(:,end))-eq_cbdc(:,2)...
    +(eq_cbdc(:,2)./(1+eq_cbdc(:,4))-eq_cbdc(:,end))*(1+i_reserve)/(1+pi_m);
output_index=Output_noCBDC/Output_noCBDC(1);
output_indexCBDC=Output_CBDC/Output_CBDC(1);
output_index=Output_noCBDC/Output_noCBDC(1);
output_indexCBDC=Output_CBDC/Output_CBDC(1);
macro_implication=figure
psi_D = eq_final_nocbdc(:,2)./(1+eq_final_nocbdc(:,4));
psi_D_R = eq_cbdc(:,2)./(1+eq_cbdc(:,4));
Delta_grid = alpha2_grid-alpha2_grid(1);
subplot(3,3,4)
hold on
h1=plot(Delta_grid/alpha3_grid(1)*100,(psi_D./psi_D(1)-1)*100,'linewidth',2);
h2=plot(Delta_grid/alpha3_grid(1)*100,(psi_D_R./psi_D_R(1)-1)*100,'--red','linewidth',2);

xlim([0,5.5]) 
xlabel('\Delta')
ylabel('Deposit')

ax=gca
ax.XAxis.FontSize = 15;
ax.YAxis.FontSize = 15;

subplot(3,3,5)
hold on
plot(Delta_grid/alpha3_grid(1)*100,(eq_final_nocbdc(:,end)./eq_final_nocbdc(1,end)-1)*100,'linewidth',2)
plot(Delta_grid/alpha3_grid(1)*100,(eq_cbdc(:,end)./eq_final_nocbdc(1,end)-1)*100,'--red','linewidth',2)
xlim([0,5.5]) 
xlabel('\Delta') 
ylabel('Loan')
ax=gca
ax.XAxis.FontSize = 15;
ax.YAxis.FontSize = 15;

subplot(3,3,6)
hold on
plot(Delta_grid/alpha3_grid(1)*100,(output_index-1)*100,'linewidth',2)
plot(Delta_grid/alpha3_grid(1)*100,(output_indexCBDC-1)*100,'--red','linewidth',2)
xlim([0,5.5]) 
xlabel('\Delta') 
ylabel('Output')
ax=gca
ax.XAxis.FontSize = 15;
ax.YAxis.FontSize = 15;

ndrate=((1+eq_final_nocbdc(:,4))*(1+pi_m)-1)*100;
ndrateR=((1+eq_cbdc(:,4))*(1+pi_m)-1)*100;
subplot(3,3,1)
hold on
plot(Delta_grid/alpha3_grid(1)*100,((1+eq_final_nocbdc(:,4))*(1+pi_m)-1)*100,'linewidth',2)
plot(Delta_grid/alpha3_grid(1)*100,((1+eq_cbdc(:,4))*(1+pi_m)-1)*100,'--red','linewidth',2)
xlim([0,5.5]) 
xlabel('\Delta') 
ylabel('Deposit Rate')
ax=gca
ax.XAxis.FontSize = 15;
ax.YAxis.FontSize = 15;

nlrate=((1+eq_final_nocbdc(:,3))*(1+pi_m)-1)*100;
nlrateR=((1+eq_cbdc(:,3))*(1+pi_m)-1)*100;
subplot(3,3,2)
hold on
plot(Delta_grid/alpha3_grid(1)*100,((1+eq_final_nocbdc(:,3))*(1+pi_m)-1)*100,'linewidth',2)
plot(Delta_grid/alpha3_grid(1)*100,((1+eq_cbdc(:,3))*(1+pi_m)-1)*100,'--red','linewidth',2)
xlim([0,5.5]) 
xlabel('\Delta') 
ylabel('Loan Rate')
ax=gca
ax.XAxis.FontSize = 15;
ax.YAxis.FontSize = 15;

subplot(3,3,3)
hold on
plot(Delta_grid/alpha3_grid(1)*100,nlrate-ndrate,'linewidth',2)
plot(Delta_grid/alpha3_grid(1)*100,nlrateR-ndrateR,'--red','linewidth',2)
xlim([0,5.5]) 
xlabel('\Delta') 
ylabel('Spread')
ax=gca
ax.XAxis.FontSize = 15;
ax.YAxis.FontSize = 15;

hL = subplot(3,3,8);
poshL = get(hL,'position');     % Getting its position

lgd = legend(hL,[h1,h2],'No CBDC','With CBDC', 'NumColumns',2, 'fontsize',15);
set(lgd,'position',poshL);      % Adjusting legend's position
axis(hL,'off');                 % Turning its axis off

saveas(macro_implication,[location,'cashless_', name ],'epsc')              % generate figure 7 in the paper
saveas(macro_implication,[location,'cashless_', name],'png')

id = find(eq_cbdc(:,end)./eq_final_nocbdc(:,end)>1,1,'first');
fraction = Delta_grid(id)/alpha3_grid(1);

disp(['A zero-interest CBDC promotes bank intermediation if ' num2str(fraction*100) '% Type 3 transactions reject cash'])

